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Abstract 

Integral functionals of Markov processes are widely used in stochastic modeling for applications in ecology, 
evolution, infectious disease epidemiology, and operations research. The integral of a stochastic process is 
often called the "cost" or "reward" accrued by the process. Many important stochastic counting models 
can be written as general birth-death processes (BDPs), which are continuous-time Markov chains on 
the non-negative integers in which only jumps to adjacent states are allowed, and there are no jumps 
down from zero. While there has been considerable progress in understanding general BDPs, work on 
integral functionals of BDPs has been limited to simple models and moments. In this paper, we show 
how to compute the distribution function and probability density of integral functionals for any general 
BDP. The method allows routine evaluation of full probability distributions for prospective modeling in 
several popular biological models and provides a framework for computing likelihoods that may be useful 
for statistical inference in a variety of applied scenarios. We provide examples of previously intractable 
integrals of BDPs from biology, queuing theory, population genetics, finance, and infectious disease 
epidemiology. In the final example, we use concepts from operations research to show how to select a 
control parameter to obtain exact probabilistic bounds on the total cost of an epidemic. 

Keywords: General birth-death process, path integral, first passage time, total cost, epidemic 



1 Introduction 



Many important real-life applications of stochastic models can be characterized as questions about path 
integrals of Markov processes. For example, the total cost of an infectious disease epidemic is the area under 



the time trajectory of the number of infected people (Jerwood 1970 Gani and Jerwood 1972 1. In ecology, 



the total time lived by a members of a population of organisms might indicate the success or failure of 



conservation efforts (Faith 1992 Crawford and Suchard 2012a I. In finance, an investor earns dividends in 



proportion to the performance of a security over time. In operations research, managers desire the quantity 
of output produced or resource consumed over time. One reason for the practical usefulness of integral 
summaries of stochastic process realizations is revealed by the units of the integral, which is often expressed 
as (objects) x (time). For example, in queuing theory, the time integral of the number of customers in a 
queue might be expressed in units of person-hours] in cancer research, the time integral of the number of 



animals alive during a trial of a radiation therapy technique might be expressed as monkey-years (Broerse 
et all 11981); traffic engineers may be interested in the number of vehicle-hours waited in models for highway 



accident delays ( jCaver 1969). 

A wide variety of discrete- valued stochastic models can be written as birth-death processes (BDPs), a 
class of continuous-time Markov chains taking values on the non-negative integers N in which only jumps 
to adjacent states are allowed, and there are no jumps down from zero (Kendall 1948). In applications, 



researchers usually treat a general BDP as a model for counting the number of hypothetical "particles" 
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(objects, organisms, species, infected people, etc.) in a system over time, where the particles can give birth 
and die. When a general BDP X{t) is in state n, a "birth" (jump to n + 1) happens with instantaneous 
rate A„ and a "death" (jump to n — 1 for n > 0) happens with instantaneous rate /i„. The study of general 
BDPs enjoyed wide interest following the groundbreaking work of Karlin and McGregor ( 1957a|b 19581, 
who derive expressions for transition probabilities, moments, and first passage times in terms of orthogonal 
polynomials. In the classical simple linear BDP, also known as the Kendall process, the per-particle birth 
and death rates are constant: A„ = nX and /i„ = (Kendall 1949). In a general BDP, A„ and /i„ are 
arbitrary functions of n, but are assumed to be time-homogeneous. In practice, the orthogonal polynomials 
and spectral measure introduced by Karlin and McGregor ( 1957a|b ) can be extremely difficult to derive for 



general BDPs (Novozhilov et al 2006 Renshaw 2011) 



In this paper, we focus on integral functionals of BDPs. Let g : N — > [0, oo) be an arbitrary positive 
function and let 5 be a set of "taboo" or prohibited states. Suppose the initial state of the BDP is X{0) = 
i G N \ S. This paper deals with the distribution of the integral functional 



g{X{t)) dt. 



where the upper limit of integration is the first passage time 

n = inf {t : X(t) G S \ X{0) = i}. 



(1) 



(2) 



Here, Wi is a functional because it maps a realization of the stochastic process g(^X{t)) to its integral. Figure 
[l] shows an example realization of a BDP and its integral Wi with S — {0}. The left-hand side shows a BDP 
beginning at X{0) — 1, and ending at X(ti) — 0. The right-hand plot shows g(^X{t)^ over the same time 
interval, and the area under the trajectory is Wi. 

The work of Karlin and McGregor ( 1957a|b[ ) provided the first theoretical tools for working with integral 



functionals of general BDPs. Puri ( 1966 1968 ) derives the characteristic function for the joint distribution 
of simple linear BDP and its integral and gives expressions for moments and limiting distributions. In a 



series of three subsequent papers, Puri further develops the theory of integrals of general stochastic processes 
and cer t ain B DPs ( |Puri[ |1971[ |1972a|b[ ). [McNeil] ( |1970[ ) gives the first rcsuhs for general BDPs, |Gani and 



McNeil ( 1971 ) derive expressions for the joint distribution of a general BDP and its integral, and Kaplan 



(1974) provides limit theorems for integrals of simple BDPs with immigration. Ball and Stefanov (2001) 
interpret BDPs as exponential families to derive generating functions for the distribution of first passage 
times and moments of stopped reward functions. More recently, researchers have found straightforward 



methods for finding moments of integrals of general BDPs using Laplace transforms ( Hernandez-Suarez and 



Castillo-Chavez 1999 Pollett and Stefanov 2003 Pollett 2003 Gani and Swift 2008). Alongside these 



advances related to integrals of BDPs, several researchers have made progress in characterizing a wide va- 



riety of interesting quantities related to BDPs in terms of continued fractions ( 
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and Guillemin (2000) introduce a combinatorial interpretation and more continued fraction expressions for 



interesting quantities related to BDPs, and Crawford and Suchard (2012b) show how to compute transi- 



tion probabilities for general BDPs using numerical inversion of continued fraction expressions for Laplace 
transforms. 

Integrals of functionals of BDPs have garnered extensive attention for their usefulness in applications as 
well. Epidemiologists have a special interest in integrals like (flj); they call this the "total cost" of an epidemic 
1970; Gani and Jerwood 



( Jerwood 
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Gani and McNeil ( 1971 1 discusses 



applications to parking lot occupancy, Gaver ( 1969 ) discusses highway traffic delays, and Pace and Khaluf 



(2012) study applications to swarm robotics. In the operations research literature, Lefevre (1981) discusses 
optimal control of epidemics that can be expressed as general BDPs; reward/cost models are often called 
"Markov decision processes" in this field. However, most analyses of integral functionals of general BDPs 
are limited to simple analytically tractable models or focused on moments of integrals like ([T]). 

In this paper, we provide the first method for evaluating the distribution of ([T]) for an arbitrary gen- 
eral BDP. We first outline the basic theory of general birth-death processes and provide expressions for the 
Laplace transform of the transition probability, first passage time distributions, and distributions of inte- 
gral functionals of general birth-death processes. Most previous work has focused on simple, analytically 
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Figure 1: Illustration of the integral of a functional of a general birth-death process (BDP). On the left, a 
BDP begins at X(0) — 1 and ends when the process reaches the absorbing state just before time t — 2. 
On the right, Wi = /J"^ g(X{t)) dt is the area under the trajectory of g{X{t)), where g : N ^ [0,oo) is an 
arbitrary positive "reward" or "cost" function. The upper limit of integration ti is the first passage time to 
zero, beginning at X{0) ~ 1. This paper is concerned with finding the probability distribution of quantities 
like Wi over all realizations of X{t). In this example, we have used a reward function g{n) that is increasing 
in n, but this need not be the case in general. 



tractable models, systems at equilibrium, or moments, which often have closed-form solutions (see, e.g. 



Hernandez- Suarez and Castillo-Chavez (19991), but often do not capture the variability in outcomes that 
might have initially motivated researchers to construct the stochastic models in the first place. In contrast, 
we establish explicit expressions for the Laplace transform of the relevant integrals using continued fractions 
and show how to evaluate these and invert them in a computationally efficient manner. The method does not 
require knowledge of the orthogonal polynomials and spectral measure introduced by |Karlin and McGregor] 
( 1957a|b" ). These developments make possible routine computation of probability distributions of integral 
functionals of general birth-death processes, and have the potential for wide use as predictive models in 
biology, medicine, and finance. The paper therefore concludes with an extensive set of concrete examples in 
which we demonstrate the distribution of certain integral functionals of previously intractable BDP models. 



2 Mathematical theory 

Let X{t) be a general BDP with transition rates and /Jk, A; = 0, 1, . . . with fiQ = 0. Then the transition 
probabilities Pij{t) — Pr(X(t) — j \ X{0) — i) satisfy the forward system of ordinary differential equations 



dPiojt) 
dt 

dt 



= ^^lP^l{t) - AoF^o(i), and 

= \j-iPi,j-i{t) + iij+iPij^i{t) — (Aj + iij)Pij{t), 



(3) 



with j > I with Pii{0) = 1 and Pij{0) = for i ^ j (Feller, 1971). The corresponding backward equations 
are 

W + ^^3P^.J+l{t) - (A, + (4) 



dP,j{t) 



dt 



2.1 Transition probabilities 



To solve for the transition probabilities, it is advantageous to work in the Laplace domain (Karlin and 



McGregor 1957b). The presentation is based on that of Crawford et al (2011). Let 



=E[e-^*] = / e-''P,,{t) dt 
Jo 



(5) 
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be the Laplace transform and Sij = 1 if i = j and zero otherwise. Laplace transforming equation ([s]) yields 

(s + Ao)/io(s) - 6io = Aii/ii(s) 
(s + \j + f-ij)ftj{s) - 6,j = Xj-ift,j-i{s) + fJ.j+ift,j+i(s) 

after straightforward rearrangement of terms. Letting i = 0, we find that 



fooi-s) = 
foj{s) 



and 



Solving for /oo(s)i we arrive at the continued fraction expression 

/ooGs) = 



s + Ao 



AoA*i 



s + Ai + ^1 



AiAf2 



(6) 



(7) 



(8) 



Letting ai — 1, aj 



s + A2 + /i2 - • • • 

-Aj_2Mj-ii ^1=5 + Aq, and bj ^ s + Xj-i + Mj-i for j > 2, (|8| becomes 



/oo(.'^) 



ai 



bi 



a2 



fli 0-2 0-3 
bi+ h2+ 63+ 



(9) 



as 



in more succinct notation. Define the fcth convergent of fm{s) as the rational function 



fli 02 



ak_ _ 



(10) 



bi+ 62+ 

With these expressions in hand, we give a fundamental result on transition probabilities for general BDPs 
Theorem 1. The Laplace transform of the transition probability Pij{t) in a general BDP is 



his) = < 




Bj{s) Bi{s)ai+2 0-1+3 
J3i+i(s)+ bi+2+ 6i+3+ 

Bi{s) Bj{s)aj+2 aj+s 



for j < i, 



for i < j. 



(11) 



The proof of this Theorem is given in Crawford and Suchard (2012b). 
2.2 Computational inversion 

With explicit expressions for the Laplace transform of the transition probability for any general BDP in 
hand, we now turn to the task of computational inversion of ( 11 ). Three advantageous facts make numerical 
inversion possible. First, there exist stable and efficient recursive schemes for evaluating continued fractions 
to an arbitrary depth; the pioneering work by Wallis (1695) and more recent developments by Lentz (19761, 



Thompson and Barnett ( 1986 ) and popularization by Press ( 2007 1 allow routine evaluation of continued 



fractions. Second, there are good a posteriori bounds for truncation of continued fractions of this type 
which allow evaluation of the error due to stopping the evaluation at a finite depth ( Craviotto et al 1993") . 
Finally, numerical Laplace inversion for well-behaved probability densities is straightforward. We appeal 
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to the methods popularized by Abate and Whitt (1992a) and Abate and Whitt (1992b). Our approach 



involves numerical evaluation of the continued fraction and numerical Laplace inversion while controlling the 
error due to truncation of the continued fraction and discretization of the inversion integral. The method is 



described in greater detail in Crawford and Suchard (2012b I, but we briefly outline the main concepts. 
Following Abate and Whitt (1992a) and Crawford and Suchard (2012b), we approximate the integral 



above by a discrete Riemann sum as follows: 



^(-l)''^Re ( 
fc=i 



A + 2fc7rV^ 
2t 



(12) 



where Re(s) is the real part of the complex variable s. The error due to discretization in (12) is 



-kA 



(13) 



k=l 



k=l 



where the second relation is follows when Pij{t) < 1, and the last when e is small (Abate and Whitt 



1992a). To obtain < 10 we set A = 7log(10). This establishes the error due to discretization of the 



Bromwich integral for Laplace inversion. However, the expressions (11 1 are infinite continued fractions which 
can only be evaluated to finite depth in practice. We therefore seek a bound on the error due to truncation 
of the continued fraction Laplace transform. Craviotto et al ( 1993 1 demonstrate the a posteriori error bound 
for convergents of continued fractions of the form (|8|): 



/oo(s) 



< 











Im 









(14) 



whenever Im(s), the imaginary part of the complex variable s, is nonzero. The expression on the right-hand 



side of ^l^ is easy to evaluate in general since one always finds /, 



way, we can control the error due to discretization in (12 1 and truncation in (14). Finally, we note that both 



^\s) before computing /oo''('S)- In this 



( 11 ) and ( 14 1 require evaluation of ratios of denominators of continued fraction convergents. Computation of 



these ratios can be highly unstable using standard methods for computing convergents; Appendix \K\ shows 
how to evaluation this ratio is a numerically stable way. 



2.3 First passage times 

Now consider the time of first arrival of a BDP X{t) into an arbitrary set S of taboo states, and suppose 
X{0) = i € N \ S. This first passage time is defined formally as 



n = inf {t : X{t) e S \ X{Q) i}. 



(15) 



To find the relationship between first passage times and the expressions for transition probabilities discussed 
above, construct a new process Y{t) identical to X{t) except that Xj = /ij = for every j e S", so every 
state in S is absorbing. Then for this modified process, with Pij{t) = Pr(F(t) ~ j \ Y{0) — i), 



(16) 



The intuitive reason for this equality is the absorbing nature of the states in S: if Y reaches a state j € S 
at any time before t, we must have Y{t) = j. Furthermore, Y cannot visit more than one state in S, so the 
absorption events are mutually exclusive and the probability of absorption is simply the sum of the individual 
absorption probabilities. Therefore the cumulative distribution function of the first passage time into S is 
given by the sum of the transition probabilities from i to every taboo state in S for the modified process 
Y(t). This observation provides the necessary link between first passage time distributions and transition 
probabilities that we will need in order to attack integrals like ([I]). 

Before we proceed with our discussion of integrals of BDPs, we must account for the possibility that the 



process first passage time may not always exist (Karlin and McGregor 1957b). In other words, the birth 
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and death rates for a general BDP may be such that the process "runs away" to infinity in finite time. 
This is also known as explosive growth. Formally, suppose the process begins at X{0) = and there are no 
absorbing states. Renshaw (2011) shows that the expected first passage time to infinity t°° is 



E(r-|X(0) = 0)=^7r,5^ 



1 



(17) 



where tti = 1 and 



n 

fc=i 



(18) 



for i > 1. When (17) diverges, the process is non-explosive, and the first passage time from to any finite 



state j is almost surely finite. When (17) is finite, the first passage time to infinity is infinite with non-zero 
probability. We return to this topic, in the context of integral functionals, in the next section. 



2.4 Integral functionals of BDPs 

Now we consider the problem of computing the distribution of ([!]). Our emphasis on first-passage times as 
the upper limit of integration in ([T]) has two benefits. First, our analyses need not be conditional on an 
arbitrary time in the future. Second, first passage times allow us to exploit powerful analytic tools that 
establish a correspondence between transition probabilities and first-passage times, allowing us to make 
analytic progress on integrals for arbitrary well-behaved processes. Our presentation follows the outline 
given by|McNeil|(|1970|). Let 



Wi{s) = E [e 



(19) 



be the Laplace transform of Wi . Recall that the destination state and waiting time U until the first departure 
from a given state are independent in a continuous-time Markov process, so the joint probability factorizes 
as follows: 



Pr(birth, [/ = u I X{0) ^ i) ^ Pr(birth | X{0) = i) x Pt{U = u \ X{Q) = i) 



Likewise, the probability of a death at time u is 

Pr(death, U \ X{Q) = i) = ^,e-(^*+^')". 



(20) 



(21) 



If the process spends time u in state i before jumping to another state, then the accumulated reward/cost is 
ug{i). Note that if X{Q) = i G S then = 0, Wi = 0, and so Wi{s) = 1. Now by an analogous conditioning 
argument for X{0) = i ^ S, we re- write (19) as 

Wi{s) = 



E 



,-s{W,+i+ug{i)) 



Pr(birth,t/ = u\X{0)^i) du 



E 



-s(W,_i+ng(i)) 



Pr(death, U \ X{0) ^i) du 



= E[e-^'^'+i] / 
Jo 

+ E[e-^'^' 
= Wi+i{s)Xi 



DO 

oo 

g-«(s£)(i) + Ai+/ti) 



(22) 
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which gives 

(s.g(z) + \i + ^ii)wi{s) = XiWi+i{s) + fj,iWi-i{s). 
Now dividing both sides of the above by g{i), we find that 

(s + A* + n*)wi{s) = A>j+i(s) + //>,_i(s) 



(23) 



(24) 



where A* = Xi/g{i) and fi* — fii/g{i). Therefore, we see that (24) is simply the backward equation for a 



modified process with birth and death rates A* and /i* for i € N. The forward equation for the cumulative 
distribution function of Wi is therefore equivalent to ^ with the modified birth and death rates. 

This remarkable fact has a simple heuristic explanation (PoUett and Stefanov 2003). Note that we can 



decompose Wi as random sum of independent and identically distributed exponential variables 



W, = Y,g{n)T^ = Y,g{n)Y,Tn 



(25) 



fc=0 



where r„ is the total time spent (sojourn time) in state n, iV„ is the number of visits to state n, and T^k is 
Exponential(A„ + /i„) waiting time during the fcth visit to state n. Multiplying each T„fe by g{n)^ we find 
that 



X] X] Exponential(A,* + /4), 



(26) 



ri=0 fc=0 



This informal derivation demonstrates that the reward/cost accumulated until the first passage time to a 
taboo state is equivalent to a first passage time under a process with modified transition rates. Pollett (2003) 



gives the conditions, analogous to those for (17), under which this modified process explodes. We note that 



differentiation of solutions of (24) yields moments of Wi, as noted by McNeil 



by Hernandez-Suarez and Castillo-Chavez (1999), Stefanov and Wang (2000), and Pollett (2003). We refer 



1970 ) and subsequently refined 



interested readers to those papers and focus here on results for the distribution of Wi, which are more useful 
in statistical and decision applications. 



To take advantage of (24 1, we modify (15) as follows. Fix S* C N and suppose X{t) is a general BDP 
with rates {A„} and with starting state X[Q) = i S N \ S*. Suppose g{n) is a positive function defined 
for all n e N. Let Y{'w) be a general BDP with rates A* = Xn/g{n) and /i* = ^n/g{n) for all n e N \ S*, and 
K = f^n = for every ne S. Then let P*j{t) = Pr(y(i) = j \ Y{0) =i). We then have 



Pr(W, < w) = ^i^*(w). 



(27) 



If instead of the cumulative distribution function of Wi we wish to have the probability density, we could 



numerically differentiate (27). However, using the properties of the Laplace transform. 



h{w) = — PiiWi < w) 
dw 



(28) 



= 5^^-i[,s/*.(s)-i^;.(0)] (w) 

= ^^-i[,s/*(s)](^) 
jes 

where flj{s) is the Laplace transform of Pij{t), ^^^[-J denotes Laplace inversion, and ^^^(0) = for all 
j G S since we have assumed i ^ S. 
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Figure 2: Probability density of the integral of a Kendall process with A = 0.1, ^ — 0.5, starting at 
X{0) = 1,2,3,4,5 until extinction. The thick black traces are exact solutions from (29 1 and the red traces 
are the numerical solutions obtained using the method outlined in section 2.4 The numerical and exact 



solutions agree to high accuracy, but the numerical method works for any general BDP and does not require 
analytic derivation of orthogonal polynomials for the process. 



3 Applications 

In this section, we discuss several practical applications of our method. The first example illustrates that our 
numerical procedure produces results that agree with known analytic expressions. We then present several 
models for which there do not exist tractable solutions. In each case, we specify the birth and death rates 
{A„} and {/in} for n = 0, 1,2, . . ., the set of taboo states S, and the reward/cost function g{n). In every 
example, we are able to generalize previous results to obtain exact distribution functions and probability 
densities for a variety of interesting integral quantities. In the final example, we present a model for the 
spread of an infectious disease and show how to set a control parameter to obtain a guaranteed probabilistic 
bound on its total cost. 



3.1 Kendall process 



In a simple linear BDP, also known as a Kendall process (Kendall 1949), A„ = nX and /i„ = nfj,. Note 
that Ao = so 5 = {0} is absorbing and corresponds to extinction in ecological and evolutionary models. 
The Kendall process is widely used in ecology as a model for reproducing and dying organisms, evolution, 
genetics, and epidemiology. To analyze integrals of the Kendall process, let g{n) = n for all n, so Wi is the 
total time lived by all particles/organisms until extinction. Under this model. A* = A for rt > 1, Aq = 0, and 
/U* = /X with /Xq = 0, which is simply the M/M/1 queue, also known as the immigration-emigration model. 



McNeil ( 1970 1 and Pollett and Stefanov] ( 2003 1 give an explicit expression for the probability density of Wi 



for /X > A (since n is almost surely finite for all i) 



dPr(T^, <w) = -e-(^+^)^(/i/A)'/2x^ (2w^) 



dw 



(29) 



where Ii{x) is the modified Bessel function of the first kind. Under certain conditions, ( |29[ ) has a limiting 
Gompertz distribution as « — >■ oo (McNeil 1970). However, using our method, we can find the distribution 
function and density of Wi without special knowledge of the functional form ( 29 ) . Figure [2] shows the 
exact probability density of Wi for i = 1,2,3,4,5 with A = 0.1 and fi = 0.5. The thick black traces are 



the exact solutions obtained from ( 29 ) and the thin red overlaid traces are the solutions computed using 
the numerical method outlined in this paper. The exact and numerically computed solutions are equal 
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to high accuracy. The benefit of our numerical method is that it is generic - it works for any BDP with 
arbitrary rates {A„} and {fJ-n}, n — 0,1, . . .. In particular, our method does not require analytic derivation 



of the orthogonal polynomials and corresponding spectral measure required by the framework of ( Karlin and 
McGregor! |1957a|b| ). 



3.2 Queues 

Many important models in queuing theory can be written as general BDPs, and the total waiting time 
(e.g. customer-hours) during a busy period is an important measure of performance in a queue. In these 
models, customers arrive into a queue (or buffer) as a Poisson process with rate A, and waiting customers 
are removed from the queue with per-customer service rate fi. In the M/M/oo queue, also known as the 
immigration-death process, there are infinitely many servers, so the arrival and service (birth and death) 
rates are A„ = A and = n/i for n > and Xq = /iq = 0, so S = {0}. Let g{x) = x and A* = A/n, Aq = 0, 
and fj,^ = iJ, for n > and — (Iq ~ 0. The distribution of the total waiting time Wi is easily found using 
our computational method. 

In the M/M/1 queue, also known as the immigration-emigration model, there is only a single server. The 
rates are A„ = A and /i„ = /i with g{x) = x for all x and S = {0}. Then the modified transition rates are 
A* = X/n and ^* = ^/n for n > and Ag = ^ 0- McNeil (1970) gives an expression for the first passage 



time of the M/M/oo queue for use in finding the distribution of a Kendall process. Analytic solution of the 
forward equations appears intractable (Daley 1969 Daley and Jacobs Jr 1969 ), and McNeil ( 1970 1 resorts 



to calculating the mean and variance. However, we can easily compute the entire probability distribution. 

Finally, we address the M/M/c queue in which there are exactly c servers. When there are n < c 
customers in the queue, the service rate is n/j,, and when there are c or more customers, the service rate 
is c/Lt for every n > c. The rates for the queue are therefore A„ — X and /i„ ~ min{n, c}/i. [Artalejo and 
Lopez- Herrero ( 2001[ ) discuss the length of the busy period and derive recursive expressions for its moments. 
The truncated nature of the death rates makes analytic solution very difficult, but poses no problem for 
our method. Figure [s] shows the person-time waited during a busy period in the M/M/oo and M/M/c 
queues, for c = 1,2,3,4,5, with customer arrival rate A = 2, mean service time = 1, and X{Q) = 5 
customers in the queue at the beginning of the time of observation. Again, we have the absorbing state 
5 = and g{n) = n. The solid line is the probability density of total waiting time in the M/M/oo queue, 
the dashed line is the M/M/5 queue, dotted M/M/4, dot-dash M/M/3, long dash M/M/2, and two-dash 
M/M/1. Note the starkly different dynamics of the waiting time distributions. We suggest that numerically 
computing densities such as these might be useful to managers in determining the number of servers in an 



M/M/c queue that achieves a certain desirable waiting time distribution. For example. Green et al 



discusses queuing approaches to optimizing staffing levels in a hospital emergency room, and Kaplan 
outlines an approach to staffing covert counter-terrorism agencies. 



(2006) 
(2012) 



3.3 Extinction of an allele under mutation and selection 

The Moran process of population genetics models the number of alleles of a certain type at a biallelic locus 



in a population of size N ( Moran , 1958 ) . Suppose the two alleles are called Ai and A2 and that n individuals 



currently have the Ai allele. When an individual dies, it is immediately replaced by a new organism whose 
allele is drawn randomly from the collection of alleles in existence: Ai with probability n/N and A2 with 
probability {N — n)/N . Suppose Ai confers a selective advantage a to carriers, and A2 confers advantage 
P < a. In the Moran process with mutation, the allele chosen by a new individual can mutate from Ai to 
A2 with probability u and vice versa with probability v. 

To derive the birth and death rates for the Moran model with mutation and selection, we outline the 
possible ways of adding or removing an allele of each type to the population. In order for a new Ai individual 
to arise (n — ^ n -I- 1), one of the N — n A2 individuals must die. Then either the new individual must draw 
its allele from one of the n Ai alleles in existence without mutation, or it must draw one of the N — n A2 
alleles with mutation to Ai. Therefore, the addition rate is 



A„ = (iV - n) [an{l - u) + I3{N - n)v] 



(30) 
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Figure 3: Probability density of the total waiting time in the M/M/oo and M/M/c queues with arrival rate 
A = 2, mean service time fi — 1, and starting with X(0) — 7 customers. The solid trace is the M/M/oo 
queue, dashed M/M/5, dotted M/M/i, dot-dash M/M/3, long dash M/M/2, and two-dash M/M/1. Note 
the differing shapes of the total waiting time distributions. Numerically computed densities such as these 
may be useful in determining the optimal number of servers an organization needs in order to achieve a total 
waiting time distribution with certain desirable properties. 



for n — 0, . . . , N with A„ = when n > N . Likewise, in order for a new A2 individual to arise (n — !■ n — 1), 
one of the n Ai individuals must die. Then either the new individual must draw its allele from one of the 
N — n A2 alleles without mutation, or it must draw one of the n Ai alleles with mutation to A2. Therefore, 
the removal rate is 

fin = n[/3{N-n){l-v)+anu], (31) 

for n = 1, . . . ,N with fiQ — Xj^ = and /i„ = when n > N. The Moran model is widely used as a 
continuous-time alternative to the Wright-Fisher process in population genetics, and expressions for tran- 



sition probabilities are known in certain simplifying cases (Karlin and McGregor 1962 Donnelly 1984 1. 
However, the model with selection and mutation is intractable using traditional methods of analysis. How- 
ever, using our technique, we can compute the full probability distribution of Wi easily and quickly. 

Suppose we are interested in the total time lived by organisms carrying the Ai allele until its eventual 
extinction, under the assumption that the Ai allele is disadvantageous (a < (3). Equivalently, we could 
also follow the total time lived by organisms carrying A2 until its fixation in the population. Let the total 
population size be iV = 100 and suppose X{0) — 50. Let the mutation rates be m = and v — 0.03, so Ai 
cannot mutate to A2, but the reverse mutation is possible; note that this makes the state absorbing. We 
therefore let S — {0} and observe that T50, the time to fixation of A2, is almost surely finite. Let a = 0.3 
and /3 — 0.5 and g{n) — n for all n. We form the modified process with A* — A„/n and /i* — i^n/Ti-- Figure|4] 
shows the distribution of the total time lived with the A\ allele until its eventual extinction in the population. 

Integrals of population genetic processes may have uses in conservation and infectious disease settings. 
If carriers of a certain allele consume a resource (e.g. food) at a different rate from non-carriers, the total 



time lived by carriers of that allele is proportional to the resource consumed ( McNeil 1970| . Alternatively, 



suppose the differential alleles are carried by a bacterial pathogen that has infected a host animal. If bacteria 
carrying one of the alleles more readily infects other hosts, then the total time lived by carriers of that allele 
is related to the likelihood of transmission. 
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Figure 4: Total time lived by organisms carrying the Ai allele until its extinction under a Moran model with 
mutation and selection for different mutation probabilities. The fitness of allele is a = 0.5, the fitness of 
allele A2 is (3 = 1, there are iV = 50 individuals in the population, and there are X{0) = 25 with the Ai 
allele at the start of observation. The traces show the probability density of organism-time lived with the Ai 
allele until its eventual extinction, for different probabilities of mutation from Ai to A2 is v, with v = 0.01 
(sohd line), v = 0.02 (dashed), v = 0.03 (dotted), v = 0.04 (dash-dotted), v = 0.05 (long dash). Note that 
even very small changes in the mutation dramatically probability alter the total lifetime distribution. 



3.4 Perpetual integral options 



The standard stochastic models for the evolution of asset prices generally have continuous sample paths 
(e.g. Weiner processes). However, there is increasing interest in discrete- valued processes in a variety of 
contexts. Perhaps the most important is the increasing availability of near-real-time financial data (Engle 
2000). For example, ( Barndorff-Nielsen et al 2010) study low-latency trading data using discrete- valued 



Levy processes. Darolles et al (2000) use a BDP to model intra-day transaction prices. Kou and Kou (2003) 



model the market capitalization of "growth stocks" as a linear BDP with immigration and emigration and 
derive a variety of useful quantities related to the equilibrium and transient characteristics of the process. 
In this section, we show how a certain type of financial option can be analyzed using our framework for 
studying integral functional of general BDPs. 

An option is a contract to buy or sell a specified asset sometime in the future (Hull 2009). An American 



option may be exercised at any time before its expiry date, and a perpetual American option can be exercised 
at any time - it does not expire. A perpetual American integral option exercised at time t entitles the buyer 



to receive payment proportional to the integral of the price of the asset up to time t ( Kramkov and Mordecki 



1994): 



g{X{t)) dt 



(32) 



Usually, the reward function has the form g{n) = an -\- bi where a, 6 > and i — X{0). The arbitrage-free 
price of an American integral option is defined as the largest expected value over all stopping rules r of the 
reward integral: 



V* = sup E 



g{X{t)) dt 



(33) 



Some authors, discussing stochastic differential equation models for asset pricing, have shown that the 
optimal time to exercise an American perpetual option is a stopping rule for the first passage time to a state 



k such that the payoff is maximized ( [Kyprianou and Pistorius[ |2003| ). That is, 



arg max 

k 




g{X{t)) dt : T,{k) = inf {t : X{t) = k \ X{0) ^ i} 



(34) 
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Figure 5: Determination of exercise price in an undiscounted arbitrage-free American perpetual integral 
option. The top panel shows the distribution of the return on a perpetual integral option, Pr(T4^i(fc) < w) for 
A: = 11, . . . , 21. The vertical line is at Wi = 10, the desired payoff. The bottom panel shows Pr(Wi(fc) < 10) 
as a function of the price k. The horizontal dashed line shows 95% probability and the vertical dashed line 
shows that k = 27 achieves the desired bound. That is, the exercise price k = 21 guarantees that the return 
will be large enough, with high probability: Pr(Wi(27) > 10) > 0.95. 



Our interest here is not the price of the option, which is given by the expectation in ( 33 ) and can be computed 
for a given k using the methods described by Hernandez-Suarez and Castillo-Chavez (19991, Stefanov and 



Wang 


(2000 


1, and 


Pollett 


(2003 



the integral for a given stopping time, we can instead seek the stopping time that guarantees a certain return 
with high probability. More formally, we wish to find the smallest exercise price k > ^(0), and therefore 
the shortest time ri(fc), that gives 

Pr(T^,(A:) < R) < a (35) 

for a desired return i? > and a small probability < a < 1, where Wi{k) is the integral with the stopping 
time Ti{k). 

Suppose X{t) is the price of an asset that moves in unit price increments (ticks) according to a BDP 
with \n = n\ -\- a and /i„ = n/i -I- (5 for n > 1 with = Q. Figure [5] shows an application. The top panel 
shows the cumulative distribution of Wi for A = 2, = 1.5, a = 0.3, and /? = 0.5, with a = and 6=1 for 
fc = 11, . . . , 21 and X{Q) = 10. The bottom panel shows Pr(iyi > 10) as a function of k. The stopping price 
= 27 is the lowest strike price to achieve Pr(Wi(fc) > 10) > 0.95. That is, Ti(27) is the shortest time we 
must wait to ensure that the return is greater than 10, with greater than 95% probability. 

We note that our approach here does not relate directly to the traditional problems of pricing of options 
or finding the optimal time to exercise. Rather, our method allows selection of an exercise strategy that 
provides a probabilistic bound on the size of the return. One drawback of our approach is that we are unable 
to provide the same result for a discounted return (e.g. e~'"^'VFi), which arises when one desires the ratio of 
return to that of a riskless investment with interest rate r. The difficulty arises because the modified process 
e*"* /q X{s) ds no longer enjoys the Markov property providing exponential waiting times. We are actively 
working on a solution for discounted processes. 
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3.5 Probabilistic control of an epidemic 



In infectious disease epidemiology, stochastic modeling can give valuable insight into both disease dynamics 



and optimal intervention strategies (Wickwire 1977 Ball 1986). To illustrate, we model the number of 



infected persons in a homogeneously mixing population as a type of general BDP. This simple model, called 
the susceptible-infected-susceptible (SIS) model, keeps track of the number of infected in a finite population 
of size N Bailey (1957). If there are currently n < N infected persons in the population, the rate of 
new infections is proportional to the product of the number of infected n a and susceptible N — n. The 
contact /transmission rate between infected and susceptible persons is A. Infected persons recover and revert 
to susceptible status with constant per-person rate fi. For a SIS process X(t), the addition and removal 
rates are 

Xn = Xn{N — n) and /i„ = -|- e) (36) 

where e is a positive control parameter related to vaccination or some other public health intervention 
strategy. Suppose the initial number of infected is X(0) = i < N and we are interested in the total cost of 
the epidemic until its eventual extinction, so S' = {0}. Let the cost of managing the epidemic per unit time 
be ae. Additionally, let the cost per infected person per unit time be 6 > 0, so the cost function becomes 
g{n) = ae + bn. Then the total cost is 



W^^ [ [ae + bX{t)] dt = aen + b [ X{t) dt 
Jq Jo 



(37) 



where Ti is the time to extinction of the epidemic. 

As in the previous examples in this section, we can easily compute the distribution of total cost Wi. 
However, this time we go a step further and seek the smallest value of the control parameter that provides 
a desired probability bound on the total cost of the epidemic. Most optimal control models seek a policy 
that minimizes the expected total cost, corresponding to the expectation of ([l]) under certain conditions on 



the intervention and cost functions (Lefevre 1981 Cai and Luo 1994 Clancy 1999 Guo and Hernandez 



Lerma 2009). But with probability distributions in hand, we can attack the problem in a subtler way. 



The availability of probability distributions for the total cost of the epidemic allows us to seek the minimal 
intervention policy that guarantees that the total cost of the epidemic is small with high probability. Let 
X{t) be the process with rates given by (36) for a certain control setting e. Then we wish to find the smallest 
e such that 

Pr (W;, < C) < 1 ~ a (38) 



where C is a desired bound on the total cost, and < a < 1 is a small probability. The left-hand side of (38) 



is simply the probability that the total cost of an epidemic is less than a certain cost, which is easily computed 
using our method. Assuming this probability is continuous and increases monotonically with e near 1 — a, it 



is straightforward to find the smallest e that satisfies (38). In our experience, these assumptions are satisfied 



for all reasonable cost functions. Additionally, this model could easily accommodate a cost function g{n) 
that incorporates the quantities n, N ~ n, and e is a much more complicated way. We emphasize that is 



possible that no e exists to satisfy (38) for certain choices of the total cost bound C and probability a. 

Figure [6] shows how to find the minimal e for a SIS process with N = 100 individuals, X{0) = 50, 
infectivity A = 0.1, recovery rate /x = 8, control cost a = 0.1, and per-infected cost b = 0.3 per unit time. 
The top traces show the cumulative distribution function of the total cost for e = 0, 0.5, 1, 1.5, 2. The vertical 
gray line shows Wi = 7, and we wish to keep the total cost less than 7 with probability 1 — a = 0.95. The 
bottom trace shows Pr(W^i < 7) as a function of e. The horizontal gray dashed line shows 0.95 probability, 
and the vertical gray dashed line shows the smallest value of e that achieves this bound. The result is 
e « 3.4. It is interesting to note that this procedure can be completely automated and does not require 
analytic insight into the SIS process beyond specification of the rates A and \i. Most traditional epidemic 



control methods (i.e. Lefevre (1981)) seek a policy to control the mean cost. The analysis presented here 
might be useful to public health policymakers whose aim, rather than controlling average cost, is to bound 
the chance of catastrophe. 
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Figure 6: Probabilistic control of a stochastic SIS epidemic. At top. the distribution of total epidemic cost 
Wi for different values of a control parameter e. The dashed gray vertical line is at w = 7, and we wish to 
keep Wi < 7 with high probability. At bottom, the probability that W, < 7 as a function of the control 
parameter e. The horizontal gray dashed line denotes 0.95, and the vertical dashed line is the smallest epsilon 
that achieves Pr(W^j < 7) > 0.95; this yields e « 3.4. In this way, we can easily find the smallest value of a 
control parameter that bounds the probability that the epidemic will exceed a certain threshold. 
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4 Discussion 



Development of analytic methods for general BDPs has proceeded down two parallel paths. The first is 
theoretical and relies on the pioneering work of Karlin and McGregor ( 1957a|b 1958 ) and later developments 
related to continued fraction expressions for Laplace transformed probabilities by [Guillemin and Pinchon| 
(1998 1999) and Flajolet and Guillemin (2000). This approach benefits from generality - most results 
apply to all general BDPs, but suffers from a lack of useful time-domain expressions. The second path is 
more practical and has focused on finding analytic time-domain expressions for probabilities and moments 



of simple BDPs (see, e.g. PoUett (2003)). Unfortunately, few authors have attempted to bridge this divide, 
though the work of Murphy and O'Donohoe (1975), Abate and Whitt (1999), and Crawford and Suchard 



(2012b) are notable exceptions. The main contribution of this work is to reconcile divergent approaches to 



analyzing integral functionals of general BDPs in the continued fraction framework and to present a unifying 
representation that is amenable to computational solution. 



5 Software 

Software implementations of the routines described in this paper are available from the author's website. 



A Appendix: stable evaluation of continued fraction expressions 



Evaluating the continued fraction expre ssions in Theorem [T] can be challenging. The algorithm of |Lentz| 
(1976) and subsequent improvements by Thompson and Barnett (1986) and Press (2007) are widely used 
to evaluate expressions hke ([7|. [Crawford and Suchard (2012b I apply this method to compute the Laplace 
transform of transition probabilities for general BDPs. However, it is also necessary to evaluate ratios of 
denominators of continued fraction convergents such as 



n+l 



Bn 



and 



Bm 
Bn 



(39) 



where Bn{s) is given by the denominator of (10). This appendix describes a stable numerical technique for 



computing these quantities, based on the Lentz algorithm. 

We begin with the following well-known fact about continued fraction convergents. 



Lemma 1. Both the numerator Ak and denominator Bk of (10) satisfy the same recurrence, due to Wallis 



(1695): 



Bk 



hAk-i 
bkBk-i 



akAk-2, and 
akBk-2, 



(40) 



with Aq = Q, Ai = ai, Bq = 1, and Bi ^ bi. 

Unfortunately, evaluation of rational approximations to continued fractions using Lemma [l] is subject to 



substantial roundoff error when the terms become small ( Press 



method to compute continued fraction convergents (Lentz 1976 
To approximate the value of foo{s), given by 



2007 ) . We briefly review the modified Lentz 
Thompson and Barnett] |1986[|Pressl|2007 1 



truncating at depth k, we write 



Ms) 

Bk{s) 



(41) 



which is the fcth rational approximant to /o,o(s). The modified Lentz method stabilizes the recurrence by 
instead computing the quantities 

Ck^^ and Dk^^ (42) 



so that 



Ak- 



J0,0 



Bk 



(43) 
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and we compute Ck and using Lemma [T] as follows: 

Cfe = 6fc + -^ and Dk= , , ^ ^ . (44) 

Returning to the issue of computing ratios of convergents, computing the numerator and denominator 
directly by (401 is appealing, but prone to catastrophic roundoff error. For concreteness, assume m < n, 
and in particular, n — m + j where j is a non-negative integer Then we must evaluate 

and — (45) 

The first of these is simply 1/D„_|_i. The second we can find from[l] Let 

Zm^, = (46) 

If j = 0, then we have 

= B„JB„, = 1 (47) 

If j = 1 then 

Zm+l — Bm+l/B,n — i/Dm (48) 

Now using [T] divide the expression for Bk+2 by B^. to form a new recurrence as follows: 

Bm+2 



Z, 



m+2 



Bm 

Bm+l B. 



"m+2 " 



arn+2 ^— (49) 



bm+2 h a„i+2 



Continuing along these lines, we have 



Zm+3 — b„i+3Z,n+2 + flm+S — (50) 

and in general, 

Zm+j = brn+jZrn+j-l dm+j Zm+j-2 ■ (51) 

This is a stable recurrence for computing Bj^+j / B^-^, as desired. 
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